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TWO-DIMENSIONAL  QUICKEST 


Solution  of  the  Depth-Averaged  Transport-Dispersion  Equation 


PART  I :  INTRODUCTION 


1.  The  mathematical  modeling  of  a  water  quality  constituent  in 
a  water  body  frequently  employs  a  spatial  averaging  in  the  vertical  di¬ 
mension  to  yield  a  depth-averaged  two-dimensional  transport-dispersion 
equation: 


3(»h) 

at 


.  a(u4>h)  .  a(v(j»h)  _  a  4.  a*h\  .  a  L  a<t>h\ 
ax  +  ay  '  a^  \  x  ax  )  +  ay  \Jy 


+  s 


where 

<j>  =  water  quality  constituent  concentration 
h  =  water  depth 
t  =  time 

u,  v  =  velocity  components  in  the  x-  and  y-directions, 
respectively 

x,  y  =  two-dimensional  Cartesian  coordinate  directions 

T  ,  T  =  dispersion  coefficients  in  the  x-  and  y-directions, 
x  ^  respectively 

S  =  source  and  sink  terms  of  the  constituent 


2.  Numerous  finite  difference  schemes  have  been  applied  to  the 
solution  of  the  transport-dispersion  equation;  however,  until  recently 
low-order  spatial  and  temporal  discretization  techniques  were  predomi¬ 
nately  used  in  practical  applications.  Typically  weighted  combinations 
of  first-order  upwind  and  second-order  central  differencing  were  em¬ 
ployed.  Unfortunately,  behavioral  errors  such  as  numerical  diffusion 
associated  with  upwind  differencing  and  the  parasitic  oscillations 
characteristic  of  central  differencing  often  rendered  these  techniques 
unsuitable  for  applications  in  transport-dispersion  models.  Conse¬ 
quently,  there  is  a  need  to  progress  to  higher  order  schemes. 

3.  The  relative  merit  of  steady-state  applications  of  spatially 
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third-order  accurate  schemes  are  well  illustrated  in  the  literature 
(Chapman  and  Kuo  1981;  Han,  Humphrey,  and  Launder  1981;  Leonard, 
Leschziner,  and  McGuirk  1978;  Leonard  1979;  Leschziner  1980;  Leschziner 
an!  Rodi  1981).  The  comparisons  of  the  spatially  third-order  accurate 
QUICK  (Quadratic  Upstream  Interpolation  for  Convective  Kinematics) 
technique  (Leonard  1979)  with  upwind  and  central  differencing  show 
QUICK  to  be  far  superior  in  steady  transport  calculations  where  prac¬ 
tical  grid  spacings  are  used.  However,  Leonard,  Leschziner,  and 
McGuirk  (1978)  and  Leschziner  (1980)  have  shown  that  QUICK  suffers  from 
a  boundedness  problem  when  applied  to  transient  problems.  Leonard 
(1979),  however,  has  presented  an  approximate  temporally  third-order 
accurate  extension  of  the  QUICK  technique  called  QUICKEST  (Quadratic 
Upstream  Interpolation  for  Convective  Kinematics  with  Estimated  Stream¬ 
ing  Terms),  which  was  specifically  designed  to  address  unidirectional 
transient  transport  problems. 

4.  Davis  and  Moore  (1982)  presented  a  two-dimensional  adapta¬ 
tion  of  the  QUICKEST  scheme.  The  authors  stated  that  the  neglect  of 
the  cross-derivative  terms  had  negligible  effects  on  the  numerical 
results.  As  a  consequence,  the  two-dimensional  QUICKEST  formulation 
of  Davis  and  Moore  (1982)  is  equivalent  to  the  superposition  of  two 
unidirectional  QUICKEST  schemes.  Results  presented  in  this  report 
will  demonstrate  that  neglect  of  the  cross-derivative  terms  can 
corrupt  the  solution. 

5.  The  purpose  of  this  report  is  to  detail  the  derivation  of  a 
two-dimensional  QUICKEST  scheme  and  compare  its  performance  to  an 
existing  third-order  accurate  Lagrangian  algorithm  of  Hinstrupt,  Ke j , 
and  Kroszynski  (1977)  referred  to  in  this  report  as  the  12-POINT 
scheme.  Test  comparisons  included  both  one-  and  two-dimensional  tran¬ 
sient  transport.  Performance  criteria  examined  included  numerical 
diffusion/amplitude,  phase,  and  mass  conservation  errors. 

6.  Part  II  of  this  report  details  the  derivation  of  the  two- 
dimensional  QUICKEST.  Part  III  summarizes  the  12-POINT  scheme  used 
in  the  comparison.  Part  IV  describes  the  one-  and  two-dimensional 
test  cases  employed  and  summarizes  the  res  Its  of  the  comparisons. 


PART  II:  DERIVATION  OF  THE  TWO-DIMENSIONAL  QUICKEST 


7 .  The  two-dimensional  version  of  QUICKEST  is  an  extension  of 
the  one-dimensional  work  of  Leonard  (1979).  The  derivation  of  QUICKEST 
is  based  on  a  conservative  control  cell  formulation  and  uses  a  spatial 
six-point  upstream  weighted  interpolation  surface  with  temporal  advec- 
tive  correction  to  obtain  third-order  approximations  to  cell  and  cell 
face  averaged  quantities. 

Control  Cell  Formulation 

8.  The  depth-averaged  transport-dispersion  equation  with  source 
and  sink  terms  neglected 

a(<t>h)  .  9(u<t»h)  .  a(v<t>h)  _  a  (r  a$hV  a  L  a*h\ 

at  ax  ay  'ax  \'x  ax /+  a^  \Jy  a^/  ™ 

is  integrated  over  a  control  cell  on  a  constant  space,  square,  compu¬ 
tational  grid  (Figure  1)  and  in  time.  For  convenience,  the  depth- 
averaged  scalar,  (J»h  ,  will  henceforth  be  designated  as  <t>  .  The  exact 
integral  formulation  is  written: 

Ay/ 2  Ax/ 2 

jf  f  ^<t>n+1  "  <J>D)  dx  dy  (2a) 

-Ay/ 2  -Ax/ 2 


At  At 

=  -Ay  J  (u^  -  uL«t»L)  dt  -  Ax  J  (vT*T  -  v^g)  dt  (2b) 

0  0 


(2c) 


in  which  the  superscripts  n  and  n+1  denote  time  levels  and  the  sub¬ 
scripts  R  ,  L  ,  T  ,  and  B  denote  right,  left,  top,  and  bottom  cell 
face  averages,  respectively. 
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Ax 


Figure  1.  Computational  grid  and  definition  of  variables.  Nodal 
values  represent  cell  averages;  cell  wall  values  represent  cell 

wall  averages 
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Cell  and  Cell  Face  Averages 


9.  A  key  feature  of  QUICKEST  is  the  use  of  a  spatial  six-point 
upstream  weighted  interpolation  surface  to  spatially  estimate  cell  and 
cell  face  averages.  To  illustrate  the  procedure,  the  computation  is 
demonstrated  for  the  right  cell  face  average  using  the  information  pre¬ 
sented  in  Figure  2.  The  u  and  v  velocity  components  are  positive. 


*-1.1-1 


V  - 


■Hi 


x 

Ax 

_y_ 

Ay 


Figure  2.  Estimation  of  right  cell 
face  average  <}>^  ,  given  u,v  >  0 

10.  Specifying  a  Gauss  backward  difference  interpolation  formula 
in  the  longitudinal  £  and  transverse  r|  direction,  a  quadratic  inter 
polation  function  can  be  written  at  x  =  iAx  ,  y  =  jAy  (Hildebrand 
1956): 


<1> 


t.n 


=  i  +  4vt  + 


m  +  p  *2 


h  *  nVn  + 


4nvev 
s  1  I  n 


n(n  +  i)  fi2 


<t>. 

J  i.j 


(3) 
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where  £  =  x/Ax  and  r|  =  y/Ay  are  local  nondimensional  coordinates, 
V  is  the  backward  difference  operator,  and  6  is  the  central  dif¬ 
ferences  operator. 

11.  The  right  cell  face  average  is  then  computed  as  follows: 


d>  = 

R 


1/2 

/ 

-1/2 


*1/2, n  dn 


if*  +  <b  ) 
2  1  i,j  i+l, jJ 


x 

8 


4>. 


6 

+  JL  <j> 

24  ‘"i.j 


12.  For  negative  transverse  velocities,  the  transverse  correction 
2 

term,  6^/24  ,  remains  the  same.  However,  for  negative  longitudinal 
velocities,  the  appropriate  cell  face  average  is  computed: 


4> 


*  i  <V 


+  4>  .  .)  - 

i+l, J 


x 

8 


<t> 

i+l,j 


6 

+ 


* 


24  i+l,j 


Approximation  of  the  Left  Hand  Side  (LHS) 

13.  The  evaluation  of  the  LHS  (Equation  2a)  consists  of  approxi¬ 
mating  the  cell  average  at  each  time  level  by  a  quadratic  interpolation 
surface,  integrating  over  the  appropriate  limits,  and  spatially  approx¬ 
imating  temporal  gradients. 

14.  Substituting  the  quadratic  interpolation  function  (3)  into 
term  (2a),  the  LHS  is  written: 

Ay/2  Ax/2  1/2  1/2 

J  J  (4>n+1  -  4>n)  dx  dy  =  AxAy  J  J  (*|fn  *  ^  H/  ^  ^ 
-  Ay/2  -Ax/2  -1/2  -  1/2 


8 


PART  IV:  TEST  SIMULATIONS  AND  RESULTS 


31.  A  systematic  program  of  numerical  experimentation  was  per¬ 
formed  in  two  phases.  The  first  phase  consisted  of  the  advection  of 
a  Gaussian  distribution  in  a  uniform  velocity  field  in  one  dimension. 
The  second  phase  consisted  of  the  solid  body  rotation  of  a  Gaussian 
distribution. 


One-Dimensional  Advection 


32.  For  spatially  constant  Courant  numbers  and  no  diffusion, 
QUICKEST  and  12-POINT  are  algebraically  equivalent  in  one  dimension. 
Therefore,  results  presented  for  the  unidirectional  advection  tests 
are  applicable  to  both  QUICKEST  and  12-POINT. 

33.  The  unidirectional  advection  of  a  Gaussian  distribution  was 

conducted  at  Courant  numbers  (uAt/Ax)  of  0.25,  0.50,  0.75,  and  1.0. 

Each  test  computation  was  run  for  200  time  steps.  In  addition,  grid 

density  effects  were  investigated  by  distributing  95  percent  of  the 

initial  distribution  over  5,  9,  or  17  grid  points,  which  corresponds 

2  2  2 

to  variances  of  (Ax)  ,  4(Ax)  ,  and  l6(Ax)  ,  respectively. 

34.  Figure  5  presents  the  results  for  an  initial  Gaussian  distri- 

2 

bution  of  variance  4(Ax)  and  varying  Courant  numbers.  Examination 
of  the  figure  reveals  negligible  phase  error;  however,  a  noticeable  but 
slight  amplitude  difference  is  seen  between  Courant  numbers  0.5  and 
0.25/0.75. 

35.  The  large  amplitude  differences  observed  between  Courant 
number  1.0,  which  reproduces  the  continuum  solution  through  exact 
point-to-point  transport,  and  Courant  numbers  0.25,  0.50,  and  0.75  are 
due  to  the  grid  density  selected.  Figure  6  is  a  comparison  of  model 
simulations  performed  at  a  Courant  number  equal  to  0.5  for  varying  grid 
density.  This  figure  clearly  demonstrates  the  importance  of  select¬ 
ing  an  appropriate  grid  density. 


The  function  on  the  right  hand  side  is  approximated  by  a  12-point, 
two-dimensional  Everett  quadratic  interpolation  function: 


4>"+!  =  (1  -  r)  (1  -  s) 

1  «  J 


♦  r(l  -  s) 


s)  [l  -  rl2-  -,  r->  S2  -  ^ V  -Sj-  «21  *•  • 
L  6  x  6  yj  1,J 

fi  -  a2  -  «21  ♦”*,  • 

L  6  X  6  yj  i+l,  j 

, .  v  r  r(2  -  r)  ,2  1  -  s2  -2|.n 

L  6  x  6  yj  1,J+1 


L  1  -  r2  *2  1  -  s2  *2 \  ^n 

\  6  6  y  i+l,j+l 


2 

Where  r  and  s  are  local  coordinates  and  6  is  the  second  central 


difference  operator. 


x  -  uAt 


-  vAt 


6  ♦.  .  =  <J>.  ,  .  -  2$.  .  +  . 

x  i,j  l-l, j  i,j  i+l, J 


30.  The  second  stage  uses  second  central  differences  to  approx¬ 
imate  the  dispersion  terms  using  the  nodal  values  updated  during  the 
first  stage. 


PART  III:  THE  12 -POINT  SCHEME 


28.  The  12-POINT  scheme  uses  grid  and  variable  definitions 
equivalent  to  those  used  in  the  QUICKEST  formulation  except  that  ve¬ 
locities  are  defined  at  nodal  points;  spatially  represents  a  cell  aver¬ 
age;  and  temporally  represents  an  average  over  the  time  interval 
(0,At) . 

29.  For  each  time  step,  two  computational  stages  are  performed:' 
a  pure  advective  transport  stage  followed  by  a  dispersion  stage.  For 
the  first  stage,  the  value  of  the  scalar  at  the  node  labeled  <j»  at 
the  n+1  time  level  in  Figure  4  is  equal  to  the  scalar  at  the  point 

•  • 


vAt 

r_— j 

uAt 

$n(x-uAt,  y-vAt) 


Figure  4.  The  12-POINT  Lagrangian  procedure 

labeled  <J>  at  n  time  level.  <t>n+1  is  located  upstream  from  4>n  a 
distance  vAt  where  v  represents  the  velocity  vector.  In  algebraic 
notation,  the  relation  is  written: 


<t>n+1(x,y)  =  <t>n(x  -  uAt,  y  -  vAt) 


19 
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The  advective  and  diffusive  contributions  through  the  top  cell  face  are 


written: 


F  =  C  |  Un  .  +  4>n  \  -  A  <Dn  .  +  624>n 

T  T  2  \  X»J  2  y  i,j  6  y  i, 


j+1/2 


-  I  .  !_  62  n  +  2  n  +  V  2 

8  °y  i,j  24  Vi. j+1/2  2  y  i,j+l/2  +  2  Vi, j+1/2 


(15b) 


1  *2  .n 

+  777  6  <P.  . 

24  x  i,j 


1 

24  x^i, j+1/2 


f  P 

RT  .  ,  .  URT  *2. 

2  Vi-1/2, j+1/2  6  Vi, j+1/2 


+  VRT  fi2  *  -  a  /I  <t>n  -  ^  fiV  -^624,n 

3  xy  i, j+1/2  T\yi,j  2  y  i, j+1/2  2  xyi, j+1/2 


where 


aT  =  V  aRT  =  fRTAt 
Ax2  Ax2 


-  «  -  fTRAt 
T  "  .  2  TR  ~  .  2 

Ay  Ay 


27.  Examination  of  expressions  15a  and  15b  reveals  that  many  of 
the  finite  difference  operators  are  centered  on  cell  faces  and  not 
nodal  points.  Shifting  these  operators  upstream  one-half  increment, 
the  flux  expressions  are  written: 

fr  *  S  |l  (♦”,/  -  “5  \*Ij  ’  '  S  ('  -  cr)  sx*lj 


2™  A<pa  +(!TR  +  4R)62  n  +VlK62  n 
2  y  i ,  j  - 1  I  2  +  6  /  y  l,j  3  °xy Pi+l/2,j-l/2 


(Right  Cell  Face) 


LHS  =  AxAy  Un+!  -  ^  . 

\  i»3 


24  6x  *1+1/2, j  “  24  6y*i+l/2,j 


+  55  6x*i-l/2,j  +  55  6yVl/2,j  <Left  Cel1  Face) 


Ct  A2* 
WT  o__4> 


25  y  i,j*l/2  -  25  Vi,jtl/2  <T°>>  Cel1  Fa"> 


Cr  2  2  J 

24  Ay*i,j-l/2  +  24  ^x*i,j-l/2 J  (Bottom  Cell  Face) 

26.  Equating  the  LHS  and  right  hand  side  (RHS),  the  flux  through 
the  right  cell  face,  given  positive  u  and  v  velocity  components,  is 
written: 


FR  =  CR  2 


i  / 4>n  +  <j,n  > 

|_2  Vi,)  i+l, j, 


fs  4  *■>  +  SV 

2  x  i,j  6  Vifl/2,j 


..  1  c2.n  1  A2An  .  aR  *2^  .  aTR  .2. 

8  x  i,j  24  Ax*i+l/2,j  2  x*i+l/2,j  2  ^y*i+l/2,j 


(15a) 


+  JL  52(j>n  -  1—  524>n 

24  y  i , j  24  y  i+l/2, j 


CTR  A  An  ,  CTR  ~2.n 

2  Vi+l/2,j-l/2  6  Ay*i+l/2, j 


.  CRCTR  ,2 
+  — 5 —  o  <P. 
3  xy  i 


n  „  (a  CR  *2,n  CTR  *2  An  ) 

i+l/2, j  '  aR\Ax*i,j  "  2  Vi+l/2,j  ‘  2  6xy*i+l/2,j / 
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23.  From  the  definition  of  the  advective  transport  equation: 


3<t>  34>  3<t»  .  u  .  _ 

3t  "  "U  3x  '  V  3y  H,0,T* 


(14) 


Substituting  Equation  14  into  13: 


.  An  13  /  w\AfA  13/  3*\  *  a.  a 

Vi.j  •  2  51  V  S/  ^  ‘  2  55  V  3yJ 


Assuming  that  u  and  v  are  approximately  constant,  the  expression 
is  written: 

A  <Dn  .  -  \  uAtAx  -  \  vAtAx  |%- 

X  1.J  2  gx2  2  3x3y 


A  *n  -  I  r  6^<ba  -  ^  (j,n 

x  i,j  2  kx  1+1/2, j  2  xy  1+1/2, j 


The  diffusive  flux  through  the  right  cell  face  is  written: 


AtAyfR  /.  _  1  r  «2An  CTR  ,.2  A 

Ax  lAx*i,j  2  CR  x  i+1/2 , j  2  xy^i+1/2 


J 


Two-Dimensional  QUICKEST  Formulation 


24.  The  QUICKEST  formulation  of  the  depth- integrated  transport 
dispersion  equation  is  written: 


.n+1  An 

v.  .  =  4>. 


i.J 


iiJ 


-F  +  F  -  F  +  F 
*R  *L  X  B 


where  Fg  ,  F^  ,  F^,  ,  and  Fg  represent  advective  and  diffusive 

contributions  through  the  right,  left,  top,  and  bottom  cell  faces, 

►  respectively. 

25.  In  Equation  9  a  caret  was  used  to  denote  the  assumption  of 
approximately  constant  velocity.  It  is  necessary  to  assign  a  cell  face 
designation  in  order  to  simplify  the  advective  transport  equation. 
Rearranging  the  LHS  (Equation  9): 


ds 


/  (I)  « -  £  /  (S?)r 


Evaluating  at  r  =  j  : 


=  (Ar  +  sA  A  )f 

1/2  r  8  0,0  1/2 


=  (A  +  sA  A  )f 
K  r  r  s  o,o 


Integrating  s  on  (0,1)  : 


[(g)  ds  =  (A  ♦  i  A  A  )  f 
J  \dt/i  n  \  r  2  r  s  /  o,c 


=  a  ♦"  .  ♦  i  fa  f”*1  -  A  *”  .) 

x  x,j  2  y  x  i,j  x  i,j/ 


where 


A  _  aII  xfl 

At,.:  <!>..,  .  -  <P.  . 
X  1»J  1+l.J  1>J 


Recognizing  that 


A  ♦?  .  -  I*  to 

X  1,J  ox 


(Sp- -(©“*- it  (£)*“■ 


Equation  12  is  written: 


A  |  |-  AtAx  =  A  •  +  I  I"  ( i^AtAx 

x  i,j  2  at  \9x/  x  i,j  2  8x  \9t/ 


ava  •yAy.vyvv-  ;•  vAv'vi 

*  •  A  *1' *  A  *A**i  ‘  -  A  «■-*  ljA«  ■  *-  •«-  -  *  -k  Vl  v--  ^  *-  v  *- 


tit) 

t  -  Atj- 


t  -  0 


i - 

f0.0 


RIGHT  CELL  FACE 


7 


-- i 

fi.o 


•  • 

f0.-1  ft.-1 

- 1 _ I _ I _ 

x  ■  (i-l)Ax  x  ■  iAx  x  -  (i+1)Ax 

Figure  3.  Definition  of  interpolation  surface  to  estimate 
the  diffusive  flux  through  the  right  cell  face 

and  the  time  interval  on  t  e  (0,1)  (Figure  3).  Specifying  Newton 
forward  difference  interpolation  formulas  in  both  r  and  s  ,  a 
quadratic  interpolation  function  can  be  written  (Hildebrand  1956) : 


f 

r,s 


1  +  rA  +  r A^  +  ^  A  +  sA  + 

r  2  r  r  s  s 


f 

0,0 


where  r 


x 

Ax 


A 


is  the  forward  difference  operator. 


22.  Note  that  following  transformation  of  coordinates 


13 


Substituting  these  definitions  into  Equation  10,  the  advective  flux 
is  written: 

4y  f  dt  =  4y4xC  [♦£  *  H  (*  “  i*  *  ’  ij  ♦  •" x  M  +  r v  ff) 


+  “l"  ^  +  2UV^  +  V  ^)\ 


"4»4*cr  p  ‘  r  Vij 


Crp  „  r  At  2  n  r  At  _ 

-IB  Av4>n  +  — —  fi  d>  +  -2 —  5^4>n 

2  AyVi+l/2,j-l/2  +  2Ax2  x  1+1/2, j  2Ay2  y  1+1/2, j 


2  2 

t  CR  -2-n  ^  CTR  «2.n  .  CRCTR  .2  .n 

6  fix*i+l/2,j  6  ^y^i+1/2, j  3  6xy*i+l/2, 


where 

«l>n  -  i  (a>n  +  d»n  \  -  I  +  —  5^<j> 

R  2  \*i,j  *i+l ,j)  8  x  i,j  24  y  1 , j 

Approximation  of  the  Diffusion  Terms 

19.  To  illustrate  the  procedure,  the  computation  is  demonstrated 
for  the  right  cell  face  with  both  u  and  v  velocity  components 
positive  (Figure  3). 

20.  The  integral 


At 

ps) 


dt  =  Ayr. 


m  r 


represents  the  diffusive  flux  through  the  right  cell  face  over  the 
time  interval  (0,At). 

21.  Construct  a  two-dimensional  interpolation  function  (x,t)-» 

(r,8)  such  that  the  right  cell  face  can  be  examined  on  x  e  (-Ax/2,  Ax/2) 


I 


for  the  right  cell  face  with  both  u  and  v  velocity  components 
positive. 

17.  The  integral 


At 

Ay  /  dt 


represents  the  advective  mass  flux  through  the  right  cell  face  during  - 
the  time  interval  (0,At).  Evaluation  of  the  integral  requires  that  u^ 
and  4>R  be  known  as  function's  of  time. 

18.  As  an  approximation  assume  that  u^  is  approximately 
constant  over  the  time  interval 


AyuR  f  4>r  dt 


Expand  the  integrand  in  a  Taylor  series  about  t  =  0 


t2  82<l> 

+  +  H.O.T. 

1  atz 


At 

r 

d> 

R 

34» 

J 

+  *  2t 

0 

o 

Integrating 


*  /a.*  ^  At2  a«t>  .  At3  a2$\  _  A  A  r  .  At  a<i>  .  At2  a2* 

AyuR  lAt*R  +  —  at  +  ~T  ^2  )  "  AyAxCR  $R  +  F  3t  +  T  ^2 


From  the  definition  of  the  advective  transport  equation 


34>  3<1>  A  r  3  V  r  ^ 

St  =  '  u  II  ’  v  3y  +  rx  ^2  ry  r2 


=  u2  8i|  ♦  2uv  A  +  ,2  A  ♦  H.O.T. 

at2  ax2  3x9y  ay2 


Recognizing  that 


<L  l  a2(A  -  Sx<i>i-H/2..i  "  6x4>i-l/2,j 


ax  W2 


Ax' 


and  similar  relations  for  »  —  (~ 

ax  V 


a  Id  <t> 


ay  \Bx2J  •  and  ay  1  » 


\dy 


Equation  7  is  written: 


LHS 


=  Axiy  |*""j  -  OV  ’  fi  [-  S  (SxVl/2,j  ‘  6xVl/2,  j  ) 

-b  (s&. 


-  624> 
,j+l/2  xi 


) 

+  * 

vv 

24 

Ax  \^C 

y  1+1/2, j 


(8) 


*  S2*-  wo  4^  -  T-  fa2*-  -4.1/0  _  G2*-  •  i/o^ 
y  i-l/2, j/  Ay  \y  i,j+l/2  y  i,j-l/2y 


Recognizing  the  Courant  numbers  Cx  =  and  and  using 


the  caret  to  denote  the  approximations  in  Equation  6,  the  LHS  is 
written: 


Q 

LHS  =  AxAy  4>n+1  -  4>n  .  -  ^  Id2* 

L  X>J  24  \  * 


-  fi2<j) 

i+1/2, j  °xv 


i-l/2, j) 


C 

JL 

24 


(6x*i.j*l/2  -  ❖i.J-lrt)  -  IS  (Sfr.,/2,j  ■  5y*i-l/2,j)  (9) 


c 

.  _z 

24 


(^y^i , j+i/2  "  ^y*i, j-1/2) 


Approximation  of  the  Advective  Terms 


16.  To  illustrate  the  procedure,  the  computation  is  demonstrated 


recognizing  that 


62  4>.  .  »  Ax 
x  x,J 


2  a2* 


62  4>.  .  w  Ay 
y  i,J 


2  3^ 

ay 


,.n+l  . n  ^  ..  3<t> 

<J>  -  4>  »  At 


the  LHS  is  written: 


15.  Differentiating  and  rearranging,  the  depth-integrated 
transport-dispersion  equation  (1)  is  written: 


where  H.O.T.  represents  the  higher  order  terms.  Assume  u  and  v 
are  approximately  constant  between  adjacent  cell  faces: 

3(u4>)  34> 

— - — —  »  u  - 

8x  3x 


3(v4>) 

By 


as  V 


B0 

3y 


(6) 


Using  the  approximations  (6)  in  Equation  5  and  substituting  into  Equa¬ 
tion  4,  the  LHS  is  written: 
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36.  The  two-dimensional  solid  body  rotation  calculations  were 

conducted  using  a  symmetric  bivariate  Gaussian  distribution  of  variance 
2 

16 (Ax)  in  a  velocity  field  of  constant  angular  frequency  equal  to 
7T/1000  rad  sec  1 .  Courant  numbers  defined  at  the  center  of  mass  of  the 
distribution  ranged  from  0.34  to  0.01  which  required  200  to  6400  time 
steps . 

37.  The  results  of  the  solid  body  rotation  calculations  are 
presented  in  Figures  7  and  8,  which  compare  the  simulation  results  of 
QUICKEST  and  12-POINT  with  the  exact  solution,  respectively.  Examina¬ 
tion  of  Figures  7  and  8  reveals  that  both  12-POINT  and  QUICKEST  exhibit 
good  phase  and  amplitude  characteristics.  However,  unlike  QUICKEST, 
whose  control  cell  formulation  ensures  mass  conservation,  12-POINT  ex¬ 
hibited  a  5  percent  mass  loss.  This  mass  conservation  error  was  sub¬ 
sequently  investigated  by  repeating  the  test  simulation  with  varying 
time  steps.  Figure  9,  a  plot  of  mass  error,  in  percent,  versus  time 

-  EXACT  SOLUTION 
O  QUICKEST 


o 


Figure  7.  Comparison  of  QUICKEST  with  the 
exact  solution  for  solid  body  rotation  of 
a  bivariate  Gaussian  distribution 


step,  At  ,  clearly  illustrates  that  the  ability  of  12-POINT  to  con¬ 
serve  mass  is  strongly  dependent  upon  the  time  step  employed.  For 
completeness,  the  dependence  of  mass  conservation  error  on  grid  density 

was  investigated  by  decreasing  the  variance  of  initial  Gaussian  distri- 
2 

bution  to  4(Ax)  .  The  resulting  difference  in  mass  conservation  error 

was  less  than  1  percent. 


Neglect  of  the  Cross-Derivative  Terms 


38.  Solid  body  rotation  experiments  were  conducted  to  investigate 
the  effects  of  neglecting  the  cross-derivative  terms  in  the  formulation 
of  the  two-dimensional  QUICKEST.  Davis  and  Moore  (1982)  stated  that  the 
neglect  of  the  cross-derivative  terms  had  a  negligible  effect  on  the 
numerical  results.  The  solid  body  rotation  calculations  were  an  exten¬ 
sion  of  the  previous  two-dimensional  test  using  all  combinations  of  ini- 

2  2 

tial  variances  4(Ax)  and  l6(Ax)  and  time  steps  2.5  and  5.0  sec. 

39.  Figure  10  presents  the  comparison  of  QUICKEST  with  and  with¬ 
out  the  cross-derivative  terms  at  time  steps  of  2.5  and  5.0  sec  with 

2  2 

initial  variances  4(Ax)  and  l6(Ax)  .  Examination  of  Figure  10 
reveals  that  neglect  of  the  cross-derivative  terms  results  in  increased 
diffusion/amplitude  errors  as  grid  density  decreases  or  time  step 
increases . 


PART  V:  CONCLUSIONS 


40.  The  purpose  of  this  study  was  to  detail  the  derivation  of  a 
two-dimensional  QUICKEST  and  to  compare  its  performance  with  12-POINT, 
as  applied  to  the  solution  of  the  advective  transport  equation.  A  sys¬ 
tematic  program  of  numerical  experimentation  was  conducted  for  both 
one-dimensional  and  two-dimensional  test  cases. 

41.  Based  on  the  results  of  the  one-dimensional  test  case,  it 
may  be  concluded  that: 

a.  QUICKEST  and  12-POINT  are  algebraically  equivalent  for 
spatially  constant  Courant  numbers. 

b.  Variations  in  the  Courant  number  result  in  small  ampli¬ 
tude  differences;  however,  simulations  performed  with 
varying  grid  densities  show  that  inadequate  grid  reso¬ 
lution  results  in  large  amplitude  errors. 

c.  Phase  errors  are  negligible. 

42.  The  two-dimensional  test  results  suggest  that: 

a.  Both  QUICKEST  and  12-POINT  possess  good  phase  and  ampli¬ 
tude  characteristics. 

b.  Unlike  QUICKEST,  which  conserved  mass  under  all  test 
conditions,  the  12-POINT  scheme  exhibits  a  mass  conser¬ 
vation  error  that  varies  directly  with  the  size  of  the 
time  step  employed. 

c.  Neglect  of  the  cross-derivative  terms  results  in 
increased  diffusion/ amplitude  errors  as  grid  density 
decreases  or  time  step  increases. 

43.  Finally,  the  work  presented  herein  is  preliminary  to  the 
development  of  a  general  purpose,  depth-integrated  water  quality  model. 
Realizing  that  the  mass  conservation  property  is  a  fundamental  require¬ 
ment,  it  is  clear  that,  of  the  two  third-order  finite  difference  tech¬ 
niques  examined,  QUICKEST  is  far  superior  for  engineering  applications 
where  practical  grid  spacing  and  time  steps  are  essential. 
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